########################################################################################################
##                                                                                                    ##
##  This is the replication file for the simulations in Tsai, Tsung-han. 2014.                        ##
##  "A Bayesian Approach to Dynamic Panel Data Models with Rarely Changing Variables."                ##
##                                                                                                    ##
########################################################################################################
########################################################################################################
##  Experiment J4: The correlation between unit effects and W is varied.                              ##
##  The correlation between unit effects and X is fixed to 0.3; T=30.                                 ##
##  J=60.                                                                                             ##
########################################################################################################
## Before running, change the working directory in line 29 to the appropriate one.
##
## Require the following files:
## (1) panel_functions.R: auxiliary functions
## (2) bayes.mlm.R: The syntax for the Bayesian multilevel model
## (3) bayes.endo2.R: The syntax for the Bayesian simultaneous equation model (BSEM)
##
## Require the following packages:
## (1) ecodist
## (2) plm
## (3) arm
## (4) R2jags

## CLEAN UP
rm(list=ls());

## SET WORKING DIRECTORY
setwd("");

## LOAD REQUIRED FUNCTIONS
source("panel_functions.R");





## LOAD PACKAGES.
library(ecodist);  # GENERATE CORRELATED SAMPLES: corgen
library(plm);
library(arm);
library(R2jags);

## SET THE SEED, WHICH TELLS US WHERE THE RANDOM PROCESS STARTS.
set.seed(04102014);

## THE NUMBER OF UNITS AND TIME PERIODS
n.state <- 60; n.time <- 30;

## TRUE VALUES OF PARAMETERS
beta <- c(0.8, 0.5, 2, -1.5, -2.5, 1.8, 3);

## GENERATE ENDOGENOUS VARIABLES
exog.var <- uncor.var(n.state=n.state, t=n.time, mu.x=c(1,0), mu.z=c(2,2), sd.x=c(3,1), sd.z=c(1,2) );
x1 <- exog.var$x1; x2 <- exog.var$x2; z1 <- exog.var$z1; z2 <- exog.var$z2;

## FINAL STORAGES
corr.out <- list()

ols.est <- list();
fe.est <- list();
re.est <- list();
fevd.est <- list();
mlr.est <- list();
mlm.est <- list();
bmlm.est <- list();
bayes.est <- list();


## THE NUMBER OF REPLICATIONS
n.iter <- 100;


## SET UP THE CORRELATION BETWEEN THE UNIT EFFECTS AND REGRESSORS: x3 OR z3 (THE OTHER IS FIXED AT 0.3)
#cor.xu <- 0.8;

r <- seq(0, 0.9, 0.1);  # EXPERIMENTS FOR DIFFERENT CORRELATION PARAMETERS

for (q in 1:length(r)) {  # LOOP FOR DIFFERENT EXPERIMENTS
cor.zu <- r[q];


## STORAGE FOR ESTIMATES OF CORRELATIONS
CORR <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

## STORAGES FOR COEFFICIENTS ESTIMATES
coef.ols.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.ols.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.fe.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.fe.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.re.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.re.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.fevd.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.fevd.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.mlr.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.mlr.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.mlm.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.mlm.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bmlm.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bmlm.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bayes.matrix <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix) <- c("ly","x1","x2","x3","z1","z2","z3");


for (g in 1:n.iter) {  # LOOP FOR REPLICATIONS


## GENERATE CORRELATED VARIABLES
endog.var <- cor.var(n.state=n.state, t=n.time, mu.delta=0, sd.delta=1, cor.xu=0.3, cor.zu=cor.zu);

x3 <- endog.var$x3; z3 <- endog.var$z3
xu.corr <- endog.var$xu.corr; zu.corr <- endog.var$zu.corr;
delta <- endog.var$delta;

## GENERATE THE OUTCOME VARIABLE AND ARRANGE THE DATASET
out.data <- sim.data(n.state=n.state, t=n.time, mu=1, phi=0.8, beta=c(0.5,2,-1.5,-2.5,1.8,3), x1=x1, x2=x2, x3=x3, z1=z1, z2=z2, z3=z3, delta=delta);


########################################
## SIMULATED DATA ANALYSIS: 
########################################
## METHOD 1: POOLED OLS
pooled.ols <- lm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data);
#summary(pooled.ols)$coef[,1:2];

coef.ols <- summary(pooled.ols)$coef[-1,1];

## MATRIX OF ESTIMATES
coef.ols.matrix[g,] <- coef.ols;


########################################
## METHOD 2: FIXED EFFECTS MODEL
fe.mod <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="within"); 
#summary(fe.mod)$coef[,1:2];

coef.fe <- summary(fe.mod)$coef[,1];

## MATRIX OF ESTIMATES
coef.fe.matrix[g,] <- coef.fe;


########################################
## METHOD 3: RANDOM EFFECTS MODEL
re.mod <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="random", random.method="amemiya"); 
#summary(re.mod)$coef[,1:2];

coef.re <- summary(re.mod)$coef[-1,1];

## MATRIX OF ESTIMATES
coef.re.matrix[g,] <- coef.re;


########################################
## METHOD 4: FIXED EFFECTS VECTOR DECOMPOSITION (fevd) MODEL
# FIRST STAGE: FIT A STANDARD FIXED EFFECTS MODEL
fevd.1 <- plm(y ~ ly+x1+x2+x3+z1+z2+z3, data=out.data, index=c("country", "year"), model="within"); 
#summary(fevd.1)$coef[,1:2];

b <- coef(fevd.1); b <- matrix(b, nrow=length(b), ncol=1, byrow=FALSE);
X.bar <- out.data[, c(12:18)];
X.bar <- as.matrix(X.bar);
e <- residuals(fevd.1);

e <- matrix(e, nrow=n.state, ncol=n.time, byrow=TRUE);
ebar <- rep(rowMeans(e), each=n.time);
ebar <- matrix(ebar, nrow=n.state*n.time, ncol=1, byrow=FALSE)

u.hat <- out.data$ybar - X.bar %*% b - ebar;

# SECOND STAGE: REGRESS THE ESTIMATED UNIT EFFECTS ON THE SLOWLY CHANGING VARIABLES
fevd.2 <- lm(u.hat ~ -1+out.data$z1+out.data$z2+out.data$z3);
#summary(fevd.2)$coef[,c(1,2)];

h <- residuals(fevd.2);

# THIRD STAGE: RERUN THE FULL MODEL WITHOUT UNIT EFFECTS BUT INCLUDE h
fevd.3 <- lm(out.data$y ~ out.data$ly + out.data$x1 + out.data$x2 + out.data$x3 + out.data$z1 + out.data$z2 + out.data$z3 + h);
#summary(fevd.3)$coef[2:8,1];

coef.fevd <- summary(fevd.3)$coef[2:8,1];

## MATRIX OF ESTIMATES
coef.fevd.matrix[g,] <- coef.fevd;


########################################
## METHOD 5: MULTILEVEL MODEL BY MLE WITHOUT WITHIN-GROUP MEANS
mlm.ind <- lmer(y ~ ly+x1+x2+x3+z1+z2+z3+(1|country), data=out.data);
#summary(mlm.ind);

coef.mlm <- coef(mlm.ind)$country[1,2:8];

## MATRIX OF ESTIMATES
coef.mlr.matrix[g,] <- as.numeric(coef.mlm);


########################################
## METHOD 6: MULTILEVEL MODELING WITH GROUP MEANS
mlm.mean <- lmer(y ~ ly+x1+x2+x3+z1+z2+z3+x3bar+z3bar+(1|country), data=out.data);
#summary(mlm.ind);

coef.mlmean <- coef(mlm.mean)$country[1,2:8];

## MATRIX OF ESTIMATES
coef.mlm.matrix[g,] <- as.numeric(coef.mlmean);


########################################
## METHOD 7: BAYESIAN MULTILEVEL MODEL WITH GROUP MEANS
N <- dim(out.data)[1];
J <- length(unique(out.data$country));
country <- out.data$country

y <- out.data$y

# COVARIATES OF Y REGRESSION
X.1 <- cbind(out.data$ly, out.data$x1, out.data$x2, out.data$x3, out.data$z1, out.data$z2, out.data$z3);
M.1 <- dim(X.1)[2];

# COVARIATES OF THE SECOND LEVEL REGRESSION
Z.1 <- cbind(1, unique(out.data$x3bar), unique(out.data$z3bar) );
L <- dim(Z.1)[2];

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "Z.1", "L");

# INITIAL VALUES
inits1 <- list(
	list( tau.y=1, B=rep(0, M.1), tau.d=0.1, G=rep(-3, L) ),
	list( tau.y=0.1, B=rep(3, M.1), tau.d=0.5, G=rep(3, L) ),
	list( tau.y=0.5, B=rep(-3, M.1), tau.d=1, G=rep(0, L) )
	)

parameters1 <- c("B");

n.chains <- length(inits1)
n.iteration <- 10000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin;  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains;  # TOTAL SAMPLES SAVED


## hand-off to JAGS
## Set the working directory



bmlm <- jags(model.file="bayes.mlm.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bmlm);

## 
mcmcburnin.mcmclist <- as.mcmc(bmlm);

	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
std <- apply(mcmc.burnin, 2, sd);
Mean <- apply(mcmc.burnin, 2, mean);
coef.bmlm <- cbind(Mean[1:7], std[1:7]);

## MATRIX OF ESTIMATES
coef.bmlm.matrix[g,] <- coef.bmlm[,1];


########################################
## METHOD 8: BAYESIAN SIMULTANEOUS EQUATION MODEL

# re-run lines 219-223

# COVARIATES OF Y REGRESSION
X.1 <- cbind(1, out.data$ly, out.data$x1, out.data$x2, out.data$x3, out.data$z1, out.data$z2, out.data$z3);
M.1 <- dim(X.1)[2];

# COVARIATES OF X3 REGRESSION
X.2 <- matrix(rep(1,N),ncol=1)
L <- dim(X.2)[2];

# COVARIATES OF W3 REGRESSION
W.1 <- matrix(rep(1,N),ncol=1)
Q <- dim(W.1)[2];

K <- 3 # THE NUMBER OF REGRESSIONS

W <- diag(K)

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W");


# INITIAL VALUES
inits1 <- list(
	list(tau.y=1, tau.x3=1, tau.z3=1, B=rep(0, M.1), G=rep(0, L), H=rep(0, Q), D=array(rep(0,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.1, tau.x3=.1, tau.z3=.1, B=rep(-3, M.1), G=rep(-3, L), H=rep(-3, Q), D=array(rep(3,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.5, tau.x3=.5, tau.z3=.5, B=rep(3, M.1), G=rep(3, L), H=rep(3, Q), D=array(rep(-3,J*K), c(J,K)), Tau.D=W )
	)

parameters1 <- c("sigma.y", "B", "sigma.d", "sigma.x", "rho.dx", "rho.dz", "sigma.z")

n.chains <- length(inits1)
n.iteration <- 10000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains  # TOTAL SAMPLES SAVED


## hand-off to JAGS
## Set the working directory



bml1 <- jags(model.file="bayes.endo2.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## 
mcmcburnin.mcmclist <- as.mcmc(bml1);








	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std <- apply(mcmc.burnin, 2, sd);
Mean <- apply(mcmc.burnin, 2, mean);
coef.bayes <- cbind(Mean[2:8], std[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix[g,] <- coef.bayes[,1];


A <- c()
A[1] <- endog.var$xu.corr; 
A[2] <- endog.var$zu.corr; 
A[3:6] <- Mean[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR[g,] <- A;


}  # END OF REPLICATIONS



########################################
##  SAVE THE RESULTS INTO THE FINAL STORAGES
corr.out[[q]] <- CORR;

ols.est[[q]] <- coef.ols.matrix;
fe.est[[q]] <- coef.fe.matrix;
re.est[[q]] <- coef.re.matrix;
fevd.est[[q]] <- coef.fevd.matrix;
mlr.est[[q]] <- coef.mlr.matrix;
mlm.est[[q]] <- coef.mlm.matrix;
bmlm.est[[q]] <- coef.bmlm.matrix;
bayes.est[[q]] <- coef.bayes.matrix;


} # END OF EACH EXPERIMENT



##
length(ols.est);
dim(ols.est[[1]]);

ols.est;
fe.est;
re.est;
fevd.est;
mlr.est;
mlm.est;
bmlm.est
bayes.est;


##  SAVE THE OUTPUT OF EXPERIMENTS
setwd("/simulationJ4/");
write.csv(corr.out, "simJ4_corr.csv", row.names=FALSE);

write.csv(ols.est, "estJ4_ols.csv", row.names=FALSE);
write.csv(fe.est, "estJ4_fe.csv", row.names=FALSE);
write.csv(re.est, "estJ4_re.csv", row.names=FALSE);
write.csv(fevd.est, "estJ4_fevd.csv", row.names=FALSE);
write.csv(mlr.est, "estJ4_mlr.csv", row.names=FALSE);
write.csv(mlm.est, "estJ4_mlm.csv", row.names=FALSE);
write.csv(bmlm.est, "estJ4_bmlm.csv", row.names=FALSE);
write.csv(bayes.est, "estJ4_bayes.csv", row.names=FALSE);




################################################################################
# REPORTING RESULTS OF SIMULATIONS
################################################################################
###################################################################################
##  Experiment J4: The correlation between unit effects and W is varied.                               
##  The correlation between unit effects and X is fixed to 0.3.                                      
###################################################################################

## CLEAN UP
rm(list=ls());

##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_ols.csv");

## TRUE VALUES OF PARAMETERS
beta <- c(0.8, 0.5, 2, -1.5, -2.5, 1.8, 3);

beta.v <- rep(beta, times=10);
beta.m <- matrix(rep(beta.v, times=100), nrow=100, ncol=70, byrow=TRUE);

##  CALCULATE AVERAGED ABSOLUTE BIAS OVER TEN EXPERIMENTS
bias.est <- abs(est.mc - beta.m);

bias.m <- matrix(rep(0, 100*7), nrow=100, ncol=7);

for (i in 1:10) {
	bias.m <- bias.m + bias.est[,(1+7*(i-1)):(7+7*(i-1))];
}

colnames(bias.m) <- c("ly","x1","x2","x3","z1","z2","z3");

bias.m <- bias.m/10;

##  CALCULATE THE MEAN OF AVERAGED ABSOLUTE BIAS ACROSS 100 REPLICATIONS
bias.avg <- apply(bias.m, 2, mean);

##  CALCULATE MONTE CARLO ERRORS (MCE) OF AVERAGED ABSOLUTE BIAS
bias.mce <- apply(bias.m, 2, sd);


##  CALCULATE RMSE OVER TEN EXPERIMENTS
bias.sq <- (est.mc - beta.m)^2;

bias.sq.m <- matrix(rep(0, 100*7), nrow=100, ncol=7);

for (i in 1:10) {
	bias.sq.m <- bias.sq.m + bias.sq[,(1+7*(i-1)):(7+7*(i-1))];
}

colnames(bias.sq.m) <- c("ly","x1","x2","x3","z1","z2","z3");

bias.sq.m <- sqrt(bias.sq.m/10);


##  CALCULATE THE MEAN OF RMSE ACROSS 100 REPLICATIONS
rmse.avg <- apply(bias.sq.m, 2, mean);

##  CALCULATE MONTE CARLO ERRORS (MCE) OF RMSE
rmse.mce <- apply(bias.sq.m, 2, sd);


##
round(rbind(bias.avg, bias.mce), 3)[,c("ly", "x3", "z3")];

round(rbind(rmse.avg, rmse.mce), 3)[,c("ly", "x3", "z3")];


# OLS AAVB
            ly    x3    z3
bias.avg 0.023 0.173 0.187
bias.mce 0.002 0.012 0.024

# OLS RMSE
rmse.avg 0.024 0.184 0.217
rmse.mce 0.002 0.011 0.026


##  FOR FIGURE 1
ols.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(ols.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(ols.bias) <- seq(0, 0.9, 0.1);

ols.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(ols.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(ols.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_fe.csv");

##  RE-RUN LINES 450-493

# FE AAVB
            ly    x3   z3
bias.avg 0.002 0.020 0.048
bias.mce 0.001 0.005 0.011

# FE RMSE
rmse.avg 0.003 0.025 0.060
rmse.mce 0.001 0.006 0.014


##  FOR FIGURE 
fe.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(fe.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fe.bias) <- seq(0, 0.9, 0.1);

fe.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(fe.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fe.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_fevd.csv");

##  RE-RUN LINES 450-493

# FEVD AAVB
            ly    x3    z3
bias.avg 0.002 0.020 0.354
bias.mce 0.001 0.005 0.026

# FEVD RMSE
rmse.avg 0.003 0.025 0.379
rmse.mce 0.001 0.006 0.023


##  FOR FIGURE 
fevd.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(fevd.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fevd.bias) <- seq(0, 0.9, 0.1);

fevd.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(fevd.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(fevd.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_re.csv");

##  RE-RUN LINES 450-493

# RE AAVB
            ly    x3    z3
bias.avg 0.002 0.020 0.083
bias.mce 0.001 0.005 0.014

# RE RMSE
rmse.avg 0.003 0.025 0.101
rmse.mce 0.001 0.006 0.016


##  FOR FIGURE 
re.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(re.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(re.bias) <- seq(0, 0.9, 0.1);

re.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(re.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(re.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_mlm.csv");

##  RE-RUN LINES 450-493

# MLM WITH GROUP MEANS, AAVB
            ly    x3    z3
bias.avg 0.002 0.020 0.045
bias.mce 0.001 0.005 0.011

# MLM WITH GROUP MEANS, RMSE
rmse.avg 0.003 0.025 0.056
rmse.mce 0.001 0.006 0.013


##  FOR FIGURE 
mlm.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(mlm.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(mlm.bias) <- seq(0, 0.9, 0.1);

mlm.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(mlm.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(mlm.rmse) <- seq(0, 0.9, 0.1);


##
est.mc <- read.csv("/replication_files/simulations/simulationJ4/estJ4_bayes.csv");

##  RE-RUN LINES 450-493

# BAYESIAN AAVB
            ly    x3    z3
bias.avg 0.002 0.020 0.045
bias.mce 0.001 0.005 0.011

# BAYESIAN RMSE
rmse.avg 0.003 0.025 0.056
rmse.mce 0.001 0.006 0.013


##  FOR FIGURE 
bayes.bias <- matrix(apply(bias.est, 2, mean), nrow=10, ncol=7, byrow=TRUE);
colnames(bayes.bias) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(bayes.bias) <- seq(0, 0.9, 0.1);

bayes.rmse <- matrix(sqrt(apply(bias.sq, 2, mean)), nrow=10, ncol=7, byrow=TRUE);
colnames(bayes.rmse) <- c("ly","x1","x2","x3","z1","z2","z3");
rownames(bayes.rmse) <- seq(0, 0.9, 0.1);



###################  MAKE PLOTS
par(mfrow=c(3,2),mar=c(0,0,2,0),oma=c(4,4,2,4));

######################## ABSOLUTE BIAS OF ly FROM xcor
corr <- seq(0, 0.9, 0.1);

plot(ols.bias[,1] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(a) Absolute Bias of LDV", xlim=c(0,1.1), ylim=c(0, 0.03), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,1], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,1], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,1], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,1], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,1], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,1], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.3, 0.97, 0.97, 0.96), c(0.02, 0.007, 0.0055, 0.004));
ml.name <- c("pooled OLS", "FE, FEVD", "RE", "MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8);


######################## RMSE OF ly FROM xcor
plot(ols.rmse[,1] ~ corr, type="n", xlab="", ylab="RMSE", main="(b) RMSE of LDV", xlim=c(0,1.1), ylim=c(0, 0.03), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,1], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,1], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,1], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,1], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,1], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,1], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.3, 0.97, 0.97, 0.96), c(0.023, 0.009, 0.007, 0.005));
ml.name <- c("pooled OLS", "FE, FEVD", "RE", "MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1);


######################## ABSOLUTE BIAS OF x3 FROM xcor
plot(ols.bias[,4] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(c) Absolute Bias of x3", xlim=c(0,1.1), ylim=c(0, 0.45), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,4], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,4], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,4], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,4], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,4], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,4], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.9, 0.97, 0.9), c(0.43, 0.2, 0.05));
ml.name <- c("pooled OLS", "RE", "FE, FEVD, MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8)


######################## RMSE OF x3 FROM xcor
plot(ols.rmse[,4] ~ corr, type="n", xlab="", ylab="RMSE", main="(d) RMSE of x3", xlim=c(0,1.1), ylim=c(0, 0.45), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,4], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,4], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,4], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,4], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,4], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,4], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.92, 0.97, 0.9), c(0.43, 0.2, 0.05));
ml.name <- c("pooled OLS", "RE", "FE, FEVD, MLM, BSEM");
text(text.loc, ml.name, cex=0.8, col="black");

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1)


######################## ABSOLUTE BIAS OF w3 FROM xcor
plot(ols.bias[,7] ~ corr, type="n", xlab="", ylab="Absolute Bias", main="(e) Absolute Bias of w3", xlim=c(0,1.1), ylim=c(0, 0.4), yaxt="n", xaxt="n");

axis(side=2, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.bias[,7], lty=1, col="black", lwd=2);
lines(corr, fe.bias[,7], lty=2, col="red", lwd=2);
lines(corr, fevd.bias[,7], lty=3, col="forestgreen", lwd=2);
lines(corr, re.bias[,7], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.bias[,7], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.bias[,7], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.97, 0.99, 0.97, 0.99), c(0.07, 0.12, 0.32, 0.1));
ml.name <- c("pooled OLS", "FE, MLM", "FEVD", "BSEM, RE");
text(text.loc, ml.name, cex=0.8, col="black");

axis(side=1, las=1, tick=TRUE, cex.axis=1);

#mtext("Absolute Bias", side=2, line=2.5, outer=F, cex=.8)


######################## RMSE OF w3 FROM xcor
plot(ols.rmse[,7] ~ corr, type="n", xlab="", ylab="RMSE", main="(f) RMSE of w3", xlim=c(0,1.1), ylim=c(0, 0.4), yaxt="n", xaxt="n");

axis(side=4, las=1, tick=TRUE, cex.axis=1);

lines(corr, ols.rmse[,7], lty=1, col="black", lwd=2);
lines(corr, fe.rmse[,7], lty=2, col="red", lwd=2);
lines(corr, fevd.rmse[,7], lty=3, col="forestgreen", lwd=2);
lines(corr, re.rmse[,7], lty=4, col="goldenrod4", lwd=2);
lines(corr, mlm.rmse[,7], lty=5, col="darkslategray", lwd=2);
lines(corr, bayes.rmse[,7], lty=6, col="darkblue", lwd=2);

text.loc <- cbind(c(0.97, 0.99, 0.97, 0.99), c(0.1, 0.15, 0.36, 0.13));
ml.name <- c("pooled OLS", "FE, MLM", "FEVD", "BSEM, RE");
text(text.loc, ml.name, cex=0.8, col="black");

axis(side=1, las=1, tick=TRUE, cex.axis=1);

#mtext("RMSE", side=4, line=2.5, outer=F, cex=1)

mtext(expression(paste("corr(x3,",delta,")")), side=1, line=2.5, outer=T, cex=1);


###################################################################################
## RESULTS OF CORRELATION
###################################################################################
## CLEAN UP
rm(list=ls());

## 
xcor_corr <- read.csv("/replication_files/simulations/simulationJ4/simJ4_corr.csv");
dim(xcor_corr);

xcor_corr[1,1:6];

avg.corr <- colMeans(xcor_corr);
mavg.corr <- matrix(avg.corr, ncol=6, byrow=TRUE);
colnames(mavg.corr) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");
rownames(mavg.corr) <- seq(0,0.9,0.1);
round(mavg.corr, 3);

    tcor.dx tcor.dz rho.dx rho.dz sigma.d sigma.y
0     0.224  -0.010  0.224 -0.015   0.998   1.002
0.1   0.223   0.091  0.225  0.092   0.989   0.999
0.2   0.224   0.175  0.225  0.175   0.994   1.002
0.3   0.222   0.259  0.222  0.249   1.012   1.003
0.4   0.221   0.347  0.226  0.325   1.001   1.001
0.5   0.224   0.445  0.220  0.419   0.997   0.998
0.6   0.221   0.533  0.218  0.522   1.006   1.003
0.7   0.221   0.637  0.220  0.617   1.000   0.999
0.8   0.223   0.752  0.221  0.722   0.980   1.003
0.9   0.223   0.869  0.218  0.845   1.014   1.001


# SINCE THE DISTRIBUTIONS ARE NOT SYMMETRIC, HIGHEST POSTERIOR DENSITY (HPD) INTERVALS ARE CALCULATED
## USE "HPDint" FUNCTION
source("panel_functions.R");
hpdint <- HPDint(xcor_corr, prob=0.95);

hpdint <- t(hpdint$HPDinterval);

corr.med <- apply(xcor_corr, 2, quantile, probs=0.5);

corr.hpd <- rbind(hpdint, corr.med);

corr.com <- rbind(corr.hpd[1,], corr.hpd[3,], corr.hpd[2,]);

#corr.com <- apply(xcor_corr, 2, quantile, probs=c(0.1,0.5,0.9));
sub <- seq(1,60,by=6);

t.cor.dx <- t(corr.com[,seq(1,60,by=6)]);
t.cor.dz <- t(corr.com[,seq(2,60,by=6)]);
cor.dx <- t(corr.com[,seq(3,60,by=6)]);
cor.dz <- t(corr.com[,seq(4,60,by=6)]);

r <- seq(0,0.9,0.1);


##  FIGURE fJ41
plot(y=t.cor.dx[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(x3,",delta,")")), col="black", lwd=4, main="", axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dx[,1], r, t.cor.dx[,3], col="grey30", lwd=3);
points(y=t.cor.dx[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dx[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dx[,1], r+0.01, cor.dx[,3], col="black", lwd=3, lty=2);
points(y=cor.dx[,2], x=r+0.01, pch="+", cex=1.5, col="black");


##  FIGURE fJ42
plot(y=t.cor.dz[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(w3,",delta,")")), col="black", lwd=4, main="", axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dz[,1], r, t.cor.dz[,3], col="grey30", lwd=3);
points(y=t.cor.dz[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dz[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dz[,1], r+0.01, cor.dz[,3], col="black", lwd=3, lty=2);
points(y=cor.dz[,2], x=r+0.01, pch="+", cex=1.5, col="black");




